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(54) Title: STABLE ADAPTIVE FILTER AND METHOD 
(57) Abstract 

Stable adaptive filter and method are disclosed. The in- 
vention solves a problem of instability associated with Fast 
Affine Projection adaptive filters caused by error accumula- 
tion in an inversion process of an auto-correlation matrix. 
The Stable FAP provides updating of the adaptive filter co- 
efficients by solving at least one system of linear equations, 
whose coefficients are the auto-correlation matrix coefficients, 
by using one of the descending iterative methods having an 
inherent stability of operation due to intrinsic feedback adjust- 
ment. The results of the solution are used to update the filter 
coefficients. The above approach is applicable for any value 
of a normalized step size ranging from zero to unity. It al- 
lows either direct determining of an updating part of the filter 
coefficients without determining an inverse auto-correlation 
matrix, or determining the inverse auto-correlation matrix by 
descending iterative methods. 
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STABLE ADAPTIVE FILTER AND METHOD 

FIELD OF THE INVENTION 

This application relates to U.S. Patent Application 
5 Serial No. 09/218,428 filed on December 22, 1998, and to 
U.S. Patent Application Serial No. 09/356,041 filed on 
July 16, 1999. The present invention relates to adaptive 
filters, and in particular, to fast affine projection 
(FAP) adaptive filters providing a stability of operation, 
10 and methods of stable FAP adaptive filtering. 

BACKGROUND OF THE INVENTION 

Adaptive filtering is a digital signal processing 
technique that has been widely used in technical areas 

15 such as, e.g., echo cancellation, noise cancellation, 
channel equalization, system identification and in 
products like, e.g., network echo cancellers, acoustic 
echo cancellers for full-duplex handsfree telephones and 
audio conference systems, active noise control, data 

20 communications systems. 

The characteristics of an adaptive filter are 
determined by its adaptation algorithm. The choice of the 
adaptation algorithm in a specific adaptive filtering 
system directly affects the performance of the system. 

25 Being simple and easily stable, the normalized least 

mean square (NLMS) adaptation algorithm, being a practical 
implementation of the least mean square (L.MS) algorithm, 
is now most widely used in the industry with a certain 
degree of success. 

30 However, because of its intrinsic weakness, tho NLMS 

algorithm converges slowly with colored training i -i-vtl s 
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like the speech, an important class of signals most 
frequently encountered in many applications such as 
telecommunications. The performance of systems 
incorporating NLMS adaptive filters very often suffers 
5 from the slow convergence nature of the algorithm. Other 
known algorithms proposed so far are either too 
complicated to implement on a commercially available low- 
cost digital signal processor (DSP) or suffer from 
numerical problems. Recently, a fast affine projection 

10 (FAP) method was proposed as described in a publication by 
Steven L. Gay and Sanjeev Tavathia (Acoustic Research 
Department, AT&T Bell Laboratories) , "The Fast Affine 
Projection Algorithm," pp. 3023 - 3026, Proceedings of the 
.International Conference on Acoustics, Speech, and Signal 

15 Processing, . May 1995, Detroit, Michigan, U.S.A. The FAP is 
a simplified version of the more complicated, and 
therefore less practical, affine projection (AP) 
algorithm. With colored train signals such as the speech, 
the FAP usually converges several times faster than the 

20 NLMS, with only a marginal increase in implementation 
complexity. 

However, a stability issue has been preventing FAP from 
being used in the industry. A prior art FAP implementation 
oscillates within a short period of time even with 

25 floating-point calculations. This results from the 
accumulation of finite precision numerical errors in a 
matrix inversion process associated with the FAP. 
Researchers have been trying to solve this problem, but no 
satisfactory answer has been found so far. A remedy 

30 proposed in the publication listed above and reinforced in 
publication by Q. G. Liu, B. Champagne, and K. C. Ho 
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(Bell-Northern Research and INRS-Telecommunications , 
Universite du Quebec) , "On the Use of a Modified Fast 
Affine Projection Algorithm in Subbands for Acoustic Echo 
Cancellation," pp. 354 - 357, Proceedings of 1996 IEEE 
5 Digital Signal Processing Workshop, Loen, Norway, 
September 1996, is to periodically re-start a new 
inversion process in parallel with the old one, and to use 
it to replace the latter so as to get rid of the 
accumulated numerical errors therein. While this can be a 

10 feasible solution for high-precision DSPs such as a 
floating-point processor, it is still not suitable for 
fixed-point DSP implementations because then the finite 
precision numerical errors would accumulate so fast that 
•the re-starting period would have to be made impractically 

15 small, not to mention the extra complexity associated with 
this part of the algorithm. 

Therefore there is a need in the industry for 
development of alternative adaptive filtering methods 
which would ensure stability of operation while providing 

20 fast convergence and reliable results. 

SUMMARY OF THE INVENTION 

It is an object of the present invention to provide an 
adaptive filter and a method of adaptive filtering which 
25 would avoid the afore-mentioned problems. 

According to one aspect of the present invention 
there is provided a method of adaptive filtering, 
comprising the steps of: 

(a) determining adaptive filter coefficients; 
30 (b) defining a normalized step size; 

(c) updating the filter coefficients, comprising: 
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4 

determining auto-correlation matrix coef f icients ' from 
a reference input signal, and 

solving at least one system of linear equations whose 
coefficients are the auto-correlation matrix coefficients, 
5 the system being solved by using a descending iterative 
method having an inherent stability of its operation, the 
results of the solution being used for updating the filter 
coefficients and the number of systems of linear equations 
to be solved being dependent on the normalized step size; 
10 (d) repeating the steps (b) and (c) required number 

of times. 

Advantageously, determining of the auto-correlation 
matrix is performed recursively. The normalized step size 
•may be chosen to be equal to any value from 0 to 1 
15 depending on the application. In the majority of 
applications, it is often set to be close to unity or 
equal to. unity. Conveniently, the normalized step size is 
within a range from about 0.9 to 1.0. Another convenient 
possibility is to set the normalized step size within a 
20 range from about 0.7 to 1.0. For the normalized step size 
close to unity, the step of solving at least one system of 
linear equations comprises solving one system of linear 
equations only. Alternatively, in some applications, e.g., 
when one needs to keep misadjustment low after 
25 convergence, it is required to set the normalized step 
size substantially less than unity, e.g. less than about 
0.7. In this situation the step of solving at least one 
system of linear equations comprises solving N systems of 
linear equations, with N being a projection order. 
30 in the embodiments of the invention, a problem of 

finding the inverse of an auto-correlation matrix which is 
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inherent for other known methods, is reduced to a problem 
of solving a system of linear equations based on the auto- 
correlation matrix. The system is solved by one of 
descending iterative methods which provide inherent 
5 stability of operation due to an intrinsic feedback 
adjustment. As a result inevitable numerical errors are 
not accumulated. In first and second embodiments of the 
invention, a steepest descent and conjugate gradient 
methods are used respectively to determine the first 
10 column of the inverse auto-correlation matrix, taking into 
account that the normalized step size is close to unity. 
In a third embodiment of the invention a steepest descent 
or conjugate gradient method is used to determine 
coefficients of the inverse auto-correlation matrix by 
15 recursively solving N systems of linear equations having 
decrementing orders. It corresponds to the case of the 
normalized step size being not close to unity. The forth 
embodiment of the invention avoids determining the inverse 
of the auto-correlation matrix. Instead, a system of 
20 linear equations is solved by using a conjugate gradient 
method resulting in a solution that can be used directly 
to determine an updating part of the filter coefficients. 
Alternatively, other known descending methods, e.g. 
steepest descent, Newton's method, PART AN , quasi-Newton's 
25 method or other known iterative descending methods may 
also be used. Conveniently, the steps of the method may be 
performed by operating with real value or complex value 
numbers . 

The method described above is suitable for a variety 
30 of applications, e.g. echo cancellation, noise 
cancellation, channel equalization, system identification 
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which are widely used in products such as network echo 
cancellers, acoustic echo cancellers for full-duplex 
handsfree telephones and audio conference systems, active 
noise control systems, data communication systems. 
5 According to another aspect of the invention there is 

provided an adaptive filter, comprising: 

a filter characterized by adaptive filter coefficients; 

means for updating the filter coefficients, including 
means for setting a normalized step size, the updating 

10 means comprising: 

a correlator for determining auto-correlation matrix 
coefficients from a reference input signal, and 

a calculator for solving at least one system of linear 
•equations whose coefficients are the auto-correlation 
15 matrix coefficients, the system being solved by using a 
descending iterative method having an inherent stability 
of -its operation, the results of the solution being used 
for updating the filter coefficients and the number of 
systems of linear equations to be solved being dependent 
20 on the normalized step size. 

Advantageously, the calculator is an iterative 
calculator. Preferably, the calculator is a steepest 
descent or a conjugate gradient calculator. Alternatively, 
it may be a calculator performing a Newton's or quasi- 
25 Newton's method, a PARTAN calculator, or another known 
iterative descending calculator providing an inherent 
stability of operation. 

Conveniently, the filter and the updating means are 
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capable of operating with real numbers. Alternatively, 
they may be capable of operating with complex numbers. 

The normalized step size may be chosen to be equal 
to any value from 0 to 1 depending on the application. In 
5 the majority of applications, the adaptive filter is often 
set with the normalized step size close to unity or equal 
to unity. Conveniently, the normalized step size is within 
a range from about 0.9 to 1.0. Another convenient 
possibility is to set the normalized step size within a 

10 range from about 0.7 to 1.0. For the normalized step size 
close to unity, the calculator provides iterative solution 
of one system of linear equations only at each time 
interval. Alternatively, in some applications, e.g., when 
one needs to keep misadjustment after convergence low, it 

15 is required to set the normalized step size substantially 
less than unity, e.g. less than about 0.7. In this 
situation the calculator provides solutions of N systems 
of linear equations, with N being a projection order. 
Conveniently, due to the symmetry of the auto-correlation 

20 matrix, determining of the inverse auto-correlation matrix 
may be performed by solving N systems of linear equations 
having decrementing orders. 

The adaptive filter as described above may be used 
for echo cancellation, noise cancellation, channel 

25 equalization, system identification or other applications 
where adaptive filtering is required. 

The adaptive filter and method described above have 
an advantage over known FAP adaptive filters by providing 
a stability of operation. The problem caused by error 

30 accumulation in matrix inversion process existing in known 
FAP filters is solved in the present invention by using 
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iterative descending methods. First, the matrix inversion 
operation is reduced to a solution of a corresponding 
system of linear equations based on the auto-correlation 
matrix. Second, the iterative descending methods, used for 
5 the solution of the above system, provide an inherent 
stability of operation due to an intrinsic feedback 
adjustment. As a result, inevitable numerical errors are 
not accumulated, thus providing stability of adaptive 
filtering . 

10 

BRIEF DESCRIPTION OF THE DRAWINGS 

The invention will now be described .in greater detail 
regarding the attached drawings in which: 

Figure 1 is a block diagram of an adaptive echo 
15 cancellation system; 

Figure 2 is a block diagram of an adaptive filter 
according to the first embodiment of the invention; 

Figure 3 is a block diagram of a steepest descent 
calculator embedded in the filter of Fig. 2; 
20 Figure 4 is a block diagram of a conjugate gradient 

calculator embedded in an adaptive filter according to a 
second embodiment of the invention; 

Figure 5 is a block diagram of an adaptive filter 
according to a third embodiment of the invention; 
25 Figure 6 is a flow-chart illustrating an operation of 

a steepest descent calculator embedded in the adaptive 
filter of Fig. 5; 

Figure 7 is a flow-chart illustrating an operation of 
a conjugate gradient calculator embedded in the adaptive 
30 filter of Fig. 5; 

Figure 8 is a block diagram of an adaptive filter 
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according to a fourth embodiment of the invention; and 

Figure 9 is a block diagram of a conjugate gradient 
calculator embedded in the adaptive filter of Fig. 8. 

5 DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS 
A. CONVENTIONS IN LINEAR ALGEBRA REPRESENTATION 

In this document, underscored letters, such, as d(n) and 
X(n) , stand for column vectors, and bold-faced ones, like 
X(n), are matrices. d(n) stands for an N-l vector consisting 
10 of the N-l upper most elements of the N vector d(n) , and 
d(n) stands for an N-l vector consisting of the N-l lower 
most elements of the N vector d(n) . A superscript "T" 
stands for the transposition of a matrix or vector. 
. B« INTRODUCTION 
15 Figure 1 presents a block diagram of an adaptive echo 

cancellation system 10 with an embedded adaptive filter 
100, the echo cancellation being chosen as an exemplary 
representation of a wide class of adaptive filtering 
applications. A digitally sampled far-end reference input 
20 signal x(n) is supplied to the adaptive filter 100 and to an 
echo path 14 producing an unwanted signal u(n), the signal 
being an echo of x(n) through the echo path 14. The echo 
path 14 can be either a long electrical path, e.g. in a 
telecommunication network, or an acoustical path, e.g. in 
25 a room. An echo canceller may be used together with a 
telecomminication network switch or a speaker phone. The 
unwanted signal u(n) is mixed up with the wanted near-end 
signal s(n) in a summer 16 to produce a response signal d(n) . 
The response signal d(n) is sent to another summer 18 
30 together with an echo estimate signal y(n) generated by the 
adaptive filter 100. The summer 18 subtracts y(n) from d(n) 
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producing an output signal e(n), to be transmitted to the 
far-end. Since the echo path is constantly changing, the 
adaptive filter must be able to continuously adapt to the 
new echo path. Therefore the goal is to produce the echo 
5 estimate signal y(n) as close to u(n) as possible, so that 
the latter is largely cancelled by the former, and e(n) 
best resembles s(n) . The output signal e(n), called the error 
signal, is then transmitted to the far-end and also sent 
to the adaptive filter 100 which uses it to adjust its 
10 coefficients. 

Note that, depending on a particular application, the 
terms "far-end" and "near-end" may need to be 
interchanged. For example, with a network echo canceller 
.in a telephone terminal, x(n) in Figure 1 is actually the 
15 near-end signal to be transmitted to the far-end, and d(n). 
in Figure 1 is the signal received from the telephone loop 
connected to the far-end. Although the terminology used 
above is based on the assumption that x(n) is the far-end 
signal and d(n) is the signal perceived at the near-end, it 
20 is done solely for convenience and does not prevent the 
invention from being applied to other adaptive filter 
applications with alternate terminology. 

The following conventions in linear algebra 
representation are used throughout the text of the present 
25 patent. Underscored letters, such as d(n) and X(n), stand 
for column vectors, and bold-faced ones, like X(n), are 
matrices. d(n) stands for an N-l vector consisting of the N~ 
1 upper most elements of the N vector d(n) , and d(n) stands 
for an N-l vector consisting of the N-l lower most elements 
30 of the N vector d(n) . A superscript "T" stands for the 
transposition of a matrix or vector. 
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1. The normalized least mean square (NIjMS) filter 

The following L-dimensional column vectors are defined 
as the reference input vector and the adaptive filter 
coefficient vector respectively, where L is the length of 
5 the adaptive filter: 



10 



X(n)s 



x(n) 
x(n-l) 

x(n-L+ 1) 



and W(n) = 



w 0 (n) 
w^n) 

w L-i( n ) 



15 



20 



25 



30 



(Equation 1) 

The part for convolution and subtraction, which 
derives the output of the adaptive echo cancellation 
system, can then be expressed as 

e(n) = d(n)-y(n) = d(n) - ^ w,(n)x(n- 1) = d(n) - X T (n)W(n) 



1 = o 



(Equation 2) 

where the superscript "T" stands for transpose of a vector 



or 



matrix. The adaptation part of the method, which 
updates the coefficient vectors based on the knowledge of 
the system behavior, is 

W(n + 1) = W(n) + 2^(n)e(n)X(n) 



H(n) = 



a 



X T (n)X(n) + 6 



(Equation 3) 

In Equation (3), |l(n) is called the adaptation step size, 
which controls the rate of change to the coefficients, a 
is a normalized step size, and 5, being a small positive 
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number, prevents |i(n) f rom going too big when there is no or 
little reference signal x(n) . 

The computations required in the NLMS filter include 
2L+2 multiply and accumulate (MAC) operations and 1 divi- 
sion per sampling interval. Details about the least mean 
square (LMS) method can be found, e.g. in classical papers 
to B. Widrow, et al . , "Adaptive Noise Cancelling: Princi- 
ples and Applications," Proceedings of the IEEE, Vol. 63, 
pp. 1692 - 1716, Dec. 1975 and B. Widrow, et al. , "Sta- 
tionary and Nonstationary Learning Characteristics of the 
LMS Adaptive Filter," Proceedings of the IEEE, Vol. 64, 
pp. 1151 - 1162, Aug. 1976. 



5 



0 



2. The Affine Projection (AP) filter 

The affine projection method is a generalization of 
the NLMS method. With N being a so-called projection 
order, we define 



d(n) = 



d(n) 
d(n-l) 

_d(n-N + 1) 



, e(n) = 



e(n) 
e(n-l) 

e(n-N+l)^ 



W(n) 



X(n) = 



x(n) x(n-l) ... x(n-N+ 1) 
x(n- 1) x(n-2) ... x(n-N) 

_x(n - L + 1 ) x(n - L) ... x(n - N - L + 2) 



(Equation 4) 

where d(n) and e(n) are N vectors and X(n) is an LxN matrix. 
Usually N is much less than L, so that X(n) having more a 
"portrait" rather than a "landscape" shape. Note that e(n) 
in Equation (4) is the a priori error vector; all its ele- 
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merits, including e(n-l), e(n-N+l), depend on W(n), as indi- 
cated in Equation (5) below. 

The convolution and subtraction part of the method is 

e(n) = d(n)-X T (n)W(n) 

(Equation 5) 

where W(n) is defined in Equation (1) . The updating part of 
the method includes the following steps 

10 W(n + 1) = W(n) + aX(n)e(n) 

R(n)£(n) = e(n) or e(n) = P(n)e(n) 
P(n) = R" ! (n) 



30 



R(n) = X T (n)X(n) + SX 

(Equation 6) 

where I is the NxN identity matrix, and a and 5 play simi- 
lar roles, as described with regards to Equation 3 . a is the 
normalized step size which may have a value from 0 to 1, 
and very often is assigned a unity value. 5 is a regular- 
ization factor that prevents- R(n), the auto-correlation 
matrix, from becoming ill-conditioned or rank-deficient, 
in which case P(n) would have too big eigenvalues causing 
instability of the method- It can be seen that an NxN 
matrix inversion operation at each sampling interval is 
needed in the AP method. 

The AP method offers a good convergence property, but 
computationally is very extensive. It needs 2LN+0(N 2 ) MACs 
at each sampling interval. For example, for N equal to 5, 
which is a reasonable choice for many practical applica- 
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tions, the AP is more than 5 times as complex as the NLMS . 

3. The Fast Affine Projection (FAP) filter 

Since the AP method is impractically expensive compu- 
5 tationally, certain simplifications have been made to 
arrive at the so-called FAP method, see, e.g. US patent 
5,428,562 to Gay. Note that here the "F", for "fast",, 
means that it saves computations, not faster convergence. 
In fact by adopting these simplifications, the performance 
10 indices, including the convergence speed, will slightly 
degrade. 

Briefly, the FAP method consists of two parts: 
(a) An approximation which is shown in Equation (7) below 
.and certain simplifications to reduce the computational 
15 load. The approximation in Equation (7) uses the scaled 
posteriori errors to replace the a priori ones in Equation 
(4) : 

e(n) 
(l-a)e(n-l) 

(1 -a) 2 e(n-2) 
(l-a^-^Cn-N+l) 

(Equation 7) 

(b) The matrix inversion operation. 

The matrix inversion may be performed by using dif- 
ferent approaches. One of them is a so-called "sliding 
windowed fast recursive least square (FRLS) " approach, 
outlined in US patent 5,428,562 to Gay, to recursively 
calculate the P(n) in Eq. 6. This results in a total 
requirement of computations to be 2L+14N MACs and 5 divi- 
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e(n)« 



e(n) 
(l-oc)e(n-l) 
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sions. In another approach, the matrix inversion lemma is 
used twice to derive P(n) at sampling interval n, see, e.g. 
Q . G. Liu, B. Champagne, and K. C. Ho (Bell-Northern 
Research and INRS-Telecommunications , Universite du 
5 Quebec) , "On the Use of a Modified Fast Af fine Projection 
Algorithm in Subbands for Acoustic Echo Cancellation" , pp. 
354 - 357, Proceedings of 1996 IEEE Digital Signal Pro- 
cessing Workshop, Loen, Norway, September 1996. It assumes 
an accurate estimate P(n-l) to start with, then derives P(n) 

10 by modifying P(n-l) based on P(n-l) and knowledge of the new 
data X(n) . The total computations needed for such a FAP 
system are 2L+3N 2 + 12N MACs and 2 divisions. Compared with 
the "sliding windowed" approach, this method offers a more 
accurate estimation for P(n) because a conventional recur- 

15 sive least square (RL.S) algorithm is used, instead of a 
fast version of it which has inevitable degradations. 

Note that, it always arrives at the most accurate and 
stable solution to solve the matrix inversion problem 
directly by using classical methods. However, these meth- 

20 ods are too expensive computationally to implement on a 
real time platform. Therefore, various alternative 
approaches with much less complexity, such as the ones 
described above, are used. The above matrix inversion 
methods have no feedback adjustment. An accurate estimate 

25 of P(n) relies heavily on an accurate starting point P(n-l). 
If P(n-l) deviates from the accurate solution, the algo- 
rithm has no way of knowing that, and will still keep 
updating it based on it and the new X(n). This means that 
errors in P(n-l), if any, will very likely accumulate and be 

30 passed on to P(n), P(n+1), P(n+2) , and so on, and therefore stay 
in the system forever. When P(n) deviates from the accurate 
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value, so will the calculated £(n), as shown in Equation 
(6) . As a result, the first expression in Equation (6) 
shows that the coefficient vector W(n) will no longer be 
updated properly. That is, W(n) can be updated in wrong 
5 directions, causing the adaptive filtering system to fail. 
A proposed remedy is to periodically re-start a new inver- 
sion process, either sliding windowed FRLS or conventional 
RLS based, in parallel with the old one, and to replace 
the old one so as to get rid of the accumulated numerical 

10 errors in the latter. While this can be a feasible solu- 
tion for high-precision DSPs such as a floating-point pro- 
cessor, it is still not suitable for fixed-point DSP 
implementations because then the finite precision numeri- 
cal errors would accumulate so fast that the . re-starting 

15 period would have to be made impractically short. 

4o Stable Fast Affine Projection Filter with a nor- 
malized step size close or equal to unity 

Usually, for maximum convergence speed, the normal - 

20 ized step size a, as indicated in Equation (6), is set to 

a" value of unity, or less than but quite close to it. This 

is the case described in the publications and the US 

patent 5,428,562 cited above. It indicates that in this 

case e(n) will have only one significant element, e(n) as the 

25 very first one. Thus, the calculation for £(n) (Eq. 6) 

reduces from the product between a matrix and a vector to 

that between a vector and a scalar, i.e. 

e(n) = e(n)P(n) 



30 



(Equation 8) 

where P(n) is the very first, i.e., left most, column of 
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the matrix P(n) . Typically, a is greater than 0.9 and less 
or equal to 1 . 0 . It is also indicated in the publication 
to Q.G. Liu cited above that, even with an a slightly less 
than that range, say about 0.7, the approximation is still 
5 acceptable. Thus, one only needs to calculate N, rather 
than all the N 2 , elements of P(n) . 

In light of the above, the problem of finding P(n), the 
inverse of the auto-correlation matrix 

10 R(n) = X T (n)X(n) + 5X 

(Equation 9) 



15 



reduces to solving a set of N linear equations 

R(n)P(n) = b , b = 



(Equation 10) 

where R(n) is symmetric and positive definite according to 
its definition Equation (9), and b is an N vector with all 

20 its elements zero except the very first, which is unity. 

Although Eq. (10) is much simpler to be solved than 
the original matrix inversion problem, it is still quite 
expensive, and especially division extensive, to do that 
with classical methods like Gaussian elimination. There- 

25 fore the obtained system of linear equations is solved by 
one of iterative descending methods which provide an 
inherent stability of operation and avoid accumulation of 
numerical errors as will be described in detail below. 
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5. Stable Fast Affine Projection Filter with general step 
size 

As mentioned above, the concept described in section 
4 above, is only suitable for applications where a rela- 
5 tively large a (the one equal to unity or less than but 
very close to unity) is needed. Although a large a is 
needed in most applications, the method of adaptive fil- 
tering wouldn't be regarded as complete without addressing 
cases with smaller normilized step sizes. For example, one 
10 way of reducing the misadjus tment (steady state output 
error) after the FAP system has converged is to use a 
small a. According to Equation (6), determining of an 
updating part of the filter coefficients may be performed 
. either by a direct solving for £(n) (second line of Eq. (6) , 
15 1st formula), or by determining an inverse auto-correla- 
tion matrix (second line of Eq.(6), second formula) with 
further calculation of e(n) . Each of the above approaches 
requires to solve N systems of linear equations based on 
the auto-correlation matrix. According to the present 
20 invention, the beneficial way to do that is to use 
descending iterative methods providing stability of opera- 
tion as will be described below. 
Co PREFERRED EMBODIMENTS OF THE INVENTION 

A method of adaptive filtering implemented in an 
25 adaptive filter 100 according to the first embodiment of 
the invention includes an iterative "steepest descent" 
technique to iteratively solve the Equation (10) . 

In general, steepest descent is a technique that 
seeks the minimum point of a certain quadratic function 
30 iteratively. At each iteration (the same as sampling 
interval in our application) , it takes three steps consec- 
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utively : 

1. to find the direction in which the parameter vec- 
tor should go. This is just the negative gradient of the 
quadratic function at the current point; 

2. to find the optimum step size for the parameter 
vector updating so that it will land at the minimum point 
along the direction dictated by the above step; and 

3 . to update the parameter vector as determined above. 
By iteratively doing the above, the steepest descent 

reaches the unique minimum of the quadratic function, 

where the gradient is zero, and continuously tracks the 

minimum if it moves. Details about the steepest descent 

method can be found, for example, in a book by David G. 

Luenberger (Stanford University) , Linear and Nonlinear 

Programming, Addison-Wesley Publishing Company, 1984. 

For an adaptive filtering application, the implied 

quadratic function is as follows 

~P T (n)R(n)P(n)-P T b 
2 

(Equation 11) 

whose gradient with respect to P(n) can be easily found as 

g = R(n)P(n)-b 

(Equation 12) 

where b is defined in Equation (10) . Note that R(n) must be 
symmetric and positive definite in order for the steepest 
descent technique to be applicable, this happens to be our 
case. Seeking the minimum, where the gradient vanishes, is 
equivalent to solving Equation (10) . The steepest descent 
is also able to track the minimum point if it moves, such 
as the case with a non-stationary input signal X(n). 
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Based on the above discussion, the stable FAP (SFAP) 
method which uses the steepest descent technique includes 
the following steps: 
Initialization : 

"1, 



W(0) = 0 , X(0) = 0 , r|(0) = 0 , R(0) = 51 , a = 1 , P(0) = 



A /5 
0 



(Equation 13) 

Updating the adaptive filter coefficients in sampling 
10 interval n including: 

recursive determining of an auto-correlation matrix: 

R(n) = R(n-l) + §(n)§ T (n)-|(n-L)§ T (n-L) 



15 



20 



(Equation 14) 
where £(n) is defined in equation (23) below, and 

determining projection coefficients by solving the 
system of linear Equations (10) using the steepest descent 
technique, the projection coefficients being the coeffi- 
cients of the inverse of the auto-correlation matrix: 



g(n) = R(n)P(n-l)- 



25 



P(n) = 



g (n)g(n) 
g f (n)R(n)g(n) 



(Equation 15) 



P(n) = P(n-l)-p(n)g(n) 



(Equation 16) 



(Equation 17) 

30 and performing an adaptive filtering for updating the fil- 
ter coefficients 
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W(n) = W(n- l) + ati N _,(n- l)X(n-N) 



y(n) = W T (n)X(n) + afj T (n- l)R(n) 



e(n) = d(n)-y(n) 



(Equation 18) 



(Equation 19) 



10 



e(n) = e(n)P(n) 



ri(n) = 



0 

fl(n-l) 



+ e(n) 



(Equation 20) 



(Equation 21) 



15 



20 



25 



30 



where £(n)is 



x(n) 
x(n-l) 

x(n-N + l) 



(Equation 22) 



(Equation 23 ) 

R(n) is the first column of R(n), R(n) is an N-l vector that 
consists of the N-l lower most elements of the N vector 
R(n), and rj(n) is an N-l vector that consists of the N-l 
upper most elements of the N vector ]](n) . 

It is important to note that feedback adjustment pro- 
vided by Equations (15), (16) and (17) does not exist in 
known prior art approaches. The prior art FAP approaches 
determine P(n) based on P(n-l) and the new incoming data X(n) 
only, without examining how well a P actually approximates 
R^Ori). Therefore inevitable numerical errors will accumu- 
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late and eventually make the system collapse. The feedback 
provided by a stable descending method, used in our inven- 
tion, uses Equation (15) to examine how well P(n-l), or the 
needed part of it, approximates R -1 (n), or its corresponding 
5 part. Then the adjustments are performed in Equations (16) 
and (17) accordingly to derive P(n), or the needed part of 
it. As just mentioned, this examination is done by evalu- 
ating g(n) in Equation (15) as the feedback error. 

The three expressions shown in Equations (15) , (16) 
10 and (17) correspond to the three steps of the steepest 
descent technique discussed above. g(n) is the gradient of 
the implied quadratic function (Equation (15)), P(n) is the 
optimum step size for parameter vector adjustment, which 
.is made in Equation (17) . As follows from Table 1, the 
15 total computational requirement of the Stable FAP method 
according to the first embodiment of the invention is 
2L±2N 2 +7N-1 MACs and 1 division.. Note, that for^the steep- 
est descent technique to work adequately for the purpose 
of adaptive filtering, the projection order N has to be 
20 chosen to assure that the ' steepest descent converges 
faster than the adaptive filter coefficients do. The 
required pre-determined value of N will depend on a par- 
ticular adaptive filtering application. 

An adaptive filter 100 according to the first embodi- 
25 ment of the invention and operating in accordance with 
the method described above is shown in Figure 2. It 
includes a filter 102 characterized by adaptive filter 
coefficients W(ri)' and means 104 for updating the coeffi- 
cients, the means being set with a normalized step size a 
30 close to its maximal value, i.e. unity. The filter 102 is 
a finite impulse response (FIR) filter which receives a 
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Table 1 



10 



15 



20 



25 



30 



Equation 


Multiply and 
accumulate 
operations 


Division 


14 


2N 




15 


N 




16 


N 2 


1 


- 17 


N 




18 


Li 




19 


L+N-l 




20. 






21 


N 




22 






Total 


2L+2N 2 4-7N-1 


1 



reference input signal x(n) and an auxiliary signal f(n) (see 
Equation (33) below), used for updating the coefficients, 
and generates a provisional echo estimate signal ' PR(n) (see 
Equation (34) below). The updating means 104 includes a 
correlator 106 for recursively determining an auto-corre- 
lation signal presented in the form of auto-correlation 
matrix coefficients R(n) based on the reference input sig- 
nal x(n), and a calculator 108 for generating projection 
coefficients P(n), the projection coefficients being part 
of the coefficients of the inverse of the auto-correlation 
matrix. The calculator 108 defines projection coefficients 
by using an iterative steepest descent method having an 
inherent stability of operation as illustrated in detail 
above. The projection coefficients are used within updat- 
ing means 104 for generation the auxiliary filter adapta- 
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15 



tion signal f(n) and an echo estimate correction signal 
EC(n) (see Equation (34) below) . The latter is used 
together with the provisional echo estimate PR(n) to pro- 
duce the echo estimate signal y(n) . 
5 A convention in Fig. 2 is the use of a thick line to 

represent the propagation of a matrix or vector signal, 
i.e., with more than one component, and the use of a thin 
line to stand for a scalar signal propagation. In Fig. 2 a 
correlator 106 determines the autocorrelation matrix R(n) 
10 in accordance with the Eg. 14 using the current and past 
x(n) samples. An "u(n) calculator" 110 calculates ]}(n) based 
on Eq. 22, and as shown in Fig. 2, n(n) is not used by the 
updating means 104 until the next sampling interval. The 
filter 102 produces the convolutional sum W T (n)X(n) . Tln^Cn-l) 
is obtained from Tl N _ x (n) by putting the latter through a 
unit delay element 111, providing a delay of one sampling 
interval, and further multiplied by the step size a in a 
Multiplier 113 . The result is used for updating the adap- 
tive filter coefficients in (Eq. 18) . ^(n-l) is dot-multi- 
plied with part of R(n) by a Dot multiplier 112, and the 
result is further multiplied by a multiplier 114 with the 
step size a to form the correction term to be added to 
W T (n)X(n) by the summer 116 to form the filter output y(n) 
(Equation (19)). The summer 18 calculates the error, or 
the output, e(n), as in Equation (20) . The scalar-vector 
multiplier 118 derives £(n) in accordance with Equation 
(21) . " 

A steepest descent calculator 108 is shown in detail 
in Figure 3. Thick lines represent the propagation of a 
matrix or vector signal, i.e., with more than one compo- 
nent, and the use of a thin line stands for a scalar sig- 
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nal propagation. In the calculator 108, the auto- 
correlation matrix R(n) and the vector P(n-l) which is a part 
of the estimated inverse of R(n-l)', are multiplied in a 
Matrix-vector multiplier 130. The vector product is fur- 
5 ther subtracted by a constant vector [10 ... 0] T in a Sum- 
mer 132 to produce the gradient vector g(n), which contains 
the feedback error information about using P(n-l) as the 
estimated inverse of R(n) . This part corresponds to Equa- 
tion (15) . The squared norm of g(n) is then found by dot- 

10 multiplying g(n) with itself in a Dot multiplier 134. It is 
used as the numerator in calculating P(n) in Equation 16. A 
Matrix-vector multiplier 13 6 finds the vector product 
between the autocorrelation matrix R(n) and the gradient 
. vector g(n). This vector product is then dot-multiplied 

15 with g(n) in another Dot multiplier 13 8 to produce the 
denominator in calculating p(n) in Equation (16) . This 
denominator is reciprocated in a Reciprocator 140, and 
then further scalar-multiplied with the aforementioned 
numerator in scalar multiplier 142 to produce (3(n). This is 

20 the only place where any division operation is performed. 
Finally, P(n) is multiplied with the gradient g(n) in a sca- 
lar-vector multiplier 144 to form the correction term to 
P(n-l). This correction term is then subtracted from P(n-l) 
in a Vector Summer 146 to derive P(n) in accordance with 

25 Equation (17). P(n-l) is obtained from P(n) by using a unit 
delay element 14 8, providing a delay of one sampling 
interval . 

Two C language prototypes implementing the steepest 
descent technique according to the first embodiment of the 
30 invention have been built. The first one is a floating 
point module, and the second one is a 16-bit fixed -point 
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DSP implementation. A floating-point module simulating the 
NLMS acoustic echo canceller design in Venture, a success- 
ful full-duplex handsfree telephone terminal product by 
Nortel Networks Corporation, and a bench mark, floating- 
5 point module that repeats a prior art FAP scheme by Q. G. 
Liu, B. Champagne, and K. C. Ho (Bell-Northern Research 
and INRS-Telecommunications , Universite du Quebec), "On 
the Use of a Modified Fast Affine Projection Algorithm in 
Subbands for Acoustic Echo Cancellation," pp. 354 - 357, 

10 Proceedings of 1996 IEEE Digital. , Signal Processing Work- 
shop, Loen, Norway, September 1996, have been also imple- 
mented for comparison purposes . The following data files 
have been prepared for processing. The source ones are 
speech files with Harvard sentences (Intermediate Refer- 

15 ence System filtered or not) sampled at 8 KHz and a white 
noise file. Out of the source files certain echo files 
have been produced by filtering the source ones with cer- 
tain measured, 1200-tap, room impulse responses. These two 
sets of files act as x(n) and d(n) respectively. The major 

20 simulation results are as follows. The bench mark prior 
art floating-point FAP scheme with L=1024 and N=5, goes 
unstable at 2 '57" (2 minutes and 57 seconds, real time, 
with 8 KHz sampling rate) with speech training, but with 
certain unhealthy signs showing up after only about .2 5 

25 seconds. These signs are in the form of improper excur- 
sions of the elements of the vector P(n), first column of 
P(n) (inverse of the matrix R(n)) . The fact that it takes 
over 2 minutes from the first appearance of unhealthy 
signs to divergence, in which period the excursions of the 

30 P(n) elements become worse and worse, shows that the coef- 
ficient updating algorithm is quite tolerant of certain 
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errors in P(n) . Once simulated random quantization noises, 
which are uniformly distributed between -0.5 bit and +0.5 
bit of a 16-bit implementation, are injected into the 
matrix inversion lemma calculation, the prior art FAP sys- 
5 tern diverges in 0.6 second. 

For comparison, within the time period of our longest 
test case (7' 40"), the portions that estimate P(n), i.e., 
Eq S . (15) -(17) of the steepest descent scheme of the 
invention with the same parameters (L=1024 and N=5), 
10 always remain stable. Furthermore, the elements in the 
vector P(n) progress as expected, without any visible 
unhealthy signs like improper excursions during the entire 
7 '40" period. The output e(n) in the steepest descent 
embodiment converges approximately at the same speed as 
15 the bench mark prior art FAP and reaches the same steady 
state echo cancellation depth as the prior art FAP and 
NLMS. The SFAP according to the first embodiment of the 
invention outperforms NLMS filter; with speech training, 
it converges in about 1 second while it takes the NLMS 
20 filter about 7 to 8 seconds to do so. 

Filters of another length L=512 have also been built 
for SFAP, the prior art FAP and NLMS. As expected, they 
converge approximately twice as fast as they do for 
L=1024. 

25 Thus, the adaptive filter and method using a steepest 

descent calculator for determining the inverse matrix 
coefficients, providing a stability of adaptive filtering, 
are provided. 

A method of adaptive filtering according to a second 
30 embodiment of the present invention uses an iterative 
"conjugate gradient" technique to iteratively solve the 
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Equation (10), the corresponding calculator being shown in 
Figure 4 . 

Conjugate gradient is a technique that also seeks the 
minimum point of a certain quadratic function iteratively. 
5 Conjugate gradient is closely related to the steepest 
descent scheme discussed above. It differs from the steep- 
est decent in that it is guaranteed to reach the minimum 
in no more than N steps, with N being the order of the 
system. That is, conjugate gradient usually converges 
10 faster than the steepest descent. At each iteration (the 
same as sampling interval in out application) , the conju- 
gate gradient takes five steps consecutively: 

1. to find the gradient of the quadratic function at 
. the current point ; 
15 2 . to find the optimum factor for adjusting the 

direction vector, along which adjustment to the parameter 
_ vector will be made; 

3. to update the direction vector as determined 

above ; 

20 4 . to find the optimum step size for the parameter 

vector updating; and 

5. to update the parameter vector as determined 

above . 

Unlike the steepest descent algorithm, which simply 
25 takes the negative gradient of the quadratic function as 
the parameter vector updating direction, conjugate gradi- 
ent modifies the negative gradient to determine an opti- 
mized direction. By iteratively doing the above, the 
scheme reaches the unique minimum of the quadratic func- 
30 tion, where the gradient is zero, in no more than N steps. 
The conjugate gradient technique also continuously tracks 
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the minimum if it moves, such as the case with non-sta- 
tionary input signal x(n) . Details about the conjugate gra- 
dient algorithm can be found, for example, in a book by 
David G. Luenberger (Stanford University), Linear and Non- 
5 linear Programming, Addison-Wesley Publishing Company, 
1984. 

For an adaptive filtering application, the implied 
quadratic function is still shown in Equation (11) , whose 
gradient with respect to P(n) is also Equation (12) . Note 

10 that R(n) must be symmetric and positive definite in order 
for the conjugate gradient technique to apply, this hap- 
pens to be our case. Seeking the minimum, where the gradi- 
ent vanishes, is equivalent to solving Equation (10) . The 
conjugate gradient is also able to track the minimum point 

15 if it moves, such as the case with non-stationary input 
signal X(n) . 



according to the second embodiment, which uses the conju- 
gate gradient technique, includes the following steps: 
20 Initialization : 



Based on the above discussion, the SFAP method 



W(0) = 0 , X(0) = 0 , tj(0) = 0 , R(0) = 51 , a = 1 , P(0) 




0 



5(0) = 0 , r srs (0) = 0 , b(0) = 0 



25 



(Equation 24) 



Updating the adaptive filter coefficients in sampling 
interval n including: 

recursive determining of an auto-correlation matrix: 
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R(n) = R(n-l) + §(n)§ T (n)-§(n-L)§ T (n-L) 



10 



(Equation 25) 
where ^(n) is defined in Equation (23) above, and 
determining projection coefficients by solving the system 
of linear Equations (10) using the conjugate technique, 
the projection coefficients being first column coeffi- 
cients of the inverse of the auto-correlation matrix: 



g(n) = R(n)P(n- 1)- 



15 



y(n) = r srs (n-l)g (n)b(n-l) 



s(n) = y(n)s(n- l)-g(n) 



(Equation 26) 



(Equation 27) 



(Equation 28) 



20 



b(n) = R(n)s(n) 



(Equation 29) 



25 



r srs( n > = 1= 

§ (n)b(n) 



p(n) = -r srs (n)g I (n)s(n) 



(Equation 30) 



30 



P(n) = P(n- l) + p(n)s(n) 



(Equation 31) 
(Equation 32) 
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and performing an adaptive filtering for updating the fil- 
ter coefficients 

W(n) = BCn-D + ailN^Cn-DXdi-N) = W(n — 1 ) + f(n)X(n - N) 

(Equation 33) 

5 

y(n) = W T (n)X(n) + afj T (n - l)R(n) = PR(n) + EC(n) 

(Equation 34) 



10 



15 



20 



25 
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e(n) = d(n)-y(n) 



e(n) = e(n)P(n) 



(Equation 35) 



(Equation 36) 



Tj(n) = 



0 

Ti(n-l) 



+ e(n) 



(Equation 37) 

where R(n) is the first column of R(n), R(n) is an N-l vector 
that consists of the N-l lower most elements of the N vec- 
tor R(n), and g(n) is an N-l vector that consists of the N-l 
upper most elements of the N vector Tj(n) . 

The five expressions shown in Equations (26), (27), 
(28), (3D and (32) respectively correspond to the five 
steps of the conjugate gradient technique discussed ear- 
lier in this section. g(n) is the gradient of the implied 
quadratic function, 7(n) is the optimum factor for updating 
the direction vector s(n) . p(n) is the optimum step size for 
parameter vector adjustment, which is made in Equation 
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(32). 

As shown in Table 2,. the total computational require- 
ment of the Stable FAP method according to the second 
embodiment of the invention is 2L+2N 2 + 9N + 1 MACs and 1 
5 division. It should be also ensured that the conjugate 
gradient converges fast enough so that the adaptive filter 
coeffients converge. 

An adaptive filter according to the second embodiment 
of the invention is similar to that of the first embodi- 

10 ment shown in Figure 2 except for the calculator 108 now 
operating in accordance with the conjugate gradient tech- 
nique and being designated by numeral 208 in Figure 4. 

The conjugate gradient calculator 208 embedded in 
.the adaptive filter of the second embodiment is shown in 

15 detail in Figure 4. Thick lines represent the propagation 
of a matrix or vector signal, i.e., with more than one 
component, and the use of a thin line stands for a scalar 
signal propagation. In the calculator 208, the autocorre- 
lation matrix R(n) and the vector P(n-l), part of the esti- 

20 mated inverse of R(n-l), are multiplied in a Matrix-vector 
Multiplier 210. The resulted vector product is subtracted 
by a constant vector [1 0 ... 0} T in a Summer 212 to pro- 
duce the gradient vector g(n) , which contains the feedback 
error information about using P(n-l) as the estimated 

25 inverse of R(n). The Matrix-vector Multiplier 210 and the 
Summer 212 implement the Equation (26) above. The gradient 
g(n) is further dot-multiplied with b(n-l), an auxiliary vec- 
tor found in the last sampling interval, in a Dot Multi- 
plier 214. The resulted scalar product is multiplied by 

30 r srs (n-l) in a Multiplier 216, to produce y(n), a factor to be 
used in adjusting s(n-l), the direction vector for adjusting 
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P(n-l). r srs (n-l) is obtained from r srs (n) by putting the latter 
through a unit delay element 218, providing a delay of one 
sampling interval. Similarly, b(n-l) is obtained from b(n) by 
using another unit delay element 220. The part of the dia- 
5 gram described in this paragraph implements Equation (27) 
shown above. With 7(n) , g(n), and s(n-l) available, s(n-l) is 
then updated into s(n) by using yet another unit delay ele- 
ment 222, with a delay of one sampling interval, scalar- 
vector Multiplier 224 and Vector Summer 226 which imple- 
10 ment operations shown in Equation (28) above. Next, the 
auxiliary vector b(n), to be used in the next sampling 
interval, is calculated as the product between R(n) and s(n) 
in another Matrix-vector Multiplier 230. This implements 
Equation (29) above. The vector b(n) is then dot-multi- 
15 plied with s(n) in yet another Dot multiplier 232, and the 
scalar product is reciprocated in a Reciprocator 234, to 
produce r srs (n) (Equation (30)). This is where the only divi- 
sion operation is. By using yet another Dot Multiplier 23 6 
and a Multiplier 238, g(n) and s(n) are dot-multiplied, and 
20 the result, being a scalar product, is multiplied with 
-r srs (n) to derive P(n), thus implementing Equation (31) 
above. Once (3(n) is available, it is multiplied with s(n) in 
another scalar-vector Multiplier 24 0 to form the correc- 
tion term to P(n-l), which is then added to P(n-l) in a Vector 
25 Summer 242 in order to derive P(n) (Equation (32) above) . 

The rest of the structure of the adaptive filter, 
employing the conjugate gradient calculator 208, is simi- 
lar to that shown in Figure 2 and described above. 
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Table 2 



10 



15 



20 



25 



30 



Equation 


Multiply and 
accumulate 
operations 


Division 


! Z 3 


2N 




Z D 


N 




27 


N+l 




28 


N 




29 


N 2 




30 


N 


1 


31 


N+l 


• 


3 2 


TvT 




33 


L 




34 


L+N-l* 




35 






36 


N 




37 






Total 


2L+2N 2 +9N+1 


1 



A C language prototype for 16-bit fixed-point DSP 
implementation of the SFAP using the conjugate gradient 
scheme has been built and studied. It has the same parame- 
ters (L=1024 and N=5) and uses same data files as the 
steepest descent prototype described above. It behaves 
very similarly to its floating-point steepest descent 
counterpart. There is no observable difference in the way 
P(n) elements progress, and they also remain stable during 
the 7 '40" longest test case period. The output e(n) in the 
conjugate gradient embodiment converges approximately at 
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the saine speed as the bench mark prior art FAP and reaches 
the same steady state echo cancellation depth as the bench 
mark prior art FAP and NLMS . The SFAP according to the 
second embodiment of the invention also ourperf ormes NLMS 
5 filter in terms of convergence speed. A conjugate gradient 
filter of another length L=512 have been also built. As 
expected, it converges twice as fast as it does for 
L=1024. 

A method of adaptive filtering according to a third 
10 embodiment of the present invention provides adaptive fil- 
tering when the normalized step size has any value from 0 
to 1 . It updates the adaptive filter coefficients by iter- 
atively solving a number of systems linear equations hav- 
ing decrementing orders to determine the inverse auto- 
15 correlation matrix in a manner described below. 

Let's prove first that, if P is the inverse of a sym- 
metric matrix R, then it is also symmetric. By definition 

RP = I , PR = I 

20 (Equation 38) 

Transposing Equation (38) we get 

P T R T = I T , R T P T = I T 

25 (Equation 39) 

respectively. Since R and I are symmetric, Equation (39) 
can be written as 

P T R = I , RP T = I 
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(Equation 40 ) 
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This means that P T is also the inverse of R. Since the 
inverse of a matrix is unique, the only possibility is 

P T = P 

5 (Equation 41) 

That is, P is symmetric. 

Based on the understanding that the inverse of a sym- 
metric matrix is also symmetric, let's consider a sampling 
interval n where we need to find an N-th order, square 
10 matrix P(n) so that 

R(n)P(n) = I 

(Equation 42) 

Equation (42) can be written in ~a scalar form 

15 X^ik(n)P kj (n) = 8 ijf Vi.je [0,N-1] 

k = 0 

( Equation "43 ) 

where r ik (n) is the element of R(n) on row i and column k, and 
p k j(n) the element of P(n) on row k and column j, and is 
defined as 



20 



25 



SiJ = { i , >f i =j 

I 0 , otherwi 



otherwise 

(Equation 44) 

We first solve the set of N linear equations defined by 
j=0 in Equation (43), for {p k0 (n), k=0, 1, N-l } , i.e. 

£^ k (n)p k0 (n) = 5 i0 , Vi e [0,N-1] 
k = o 
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(Equation 45) 

Equation (45) coincides with Equation (10) derived earlier 
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and applied to the first and second embodiments of the 
invention . 



25 



R(n)P(n) = 



(Equation 46 ) 

The right hand side of Equation (45) or Equation (46) 
tells that P(n) is the left-most column of P(n) and, based 
on Equation (41), P (n) is also the upper-most row of P(n) . 

10 According to the first and second embodiments of the 
invention discussed above, this -part will cost "2N 2 + 3N" 
MACs and 1 division with steepest descent or w 2N 2 + 5N+2" 
MACs and 1 division with conjugate gradient. 

Having dealt with the j=0 case, we now start solving 

j 5 the set of N linear equations defined by j=l in Equation 
(43), for {p kl (n),k=0, 1,...,N-1}, i.e. 

N ^r ik (n)p kI (n) = 8 n , Vie [O.N-l] 

k = 0 

20 (Equation 47) 

Because P(n) is symmetric so that Poi(n) equals pjo(n), Equa- 
tion (47) can be re-arranged to become 



N X r ik( n )Picl( n ) = 5 ii- r io( n )Pio( n )> [O.N-l] 
k = 1 



(Equation 48) 

with still N equations but only N-l instead of N unknowns, 
i.e., {p k i(n), k=l, 2, N-l } , to solve. In general, these N-l 
30 unknowns can be uniquely determined by only N-l equations. 
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Thus, the equation in Equation (48) with i=0 can be omitted 
so that it becomes 

N X r ik( n )Pki( n ) = S n - r i0 (n)p 10 (n) , Vie [1,N- 1] 
k = 1 

5 

(Equation 49) 

Equation (49) has the same format as Equation (45) except 
that the order is reduced by one. Equation (49) can also 
be solved by using either of the two approaches presented 
10 above, costing "2 (N-l) 2 +4 (N-l ) MACs and 1 division with 
steepest descent" or "2 (N-l ) 2 + 6 (N-l ) +2 MACs and 1 division 
with conjugate gradient," where the added "(N-l)" in each 
of the two expressions accounts for the extra computations 
needed to calculate the right hand side of Equation (49). 
15 By repeating the above recursion steps, with the 

order of the problem decrementing by one each step, we can 
completely .solve the lower triangle of P(n). Since P(n) is 
symmetric, this is equivalent to solving the entire P(n) . A 
formula for this entire process can be derived from Equa- 
20 tion (43) and the concept described above, as follows: 
For j = 0, 1, ...,N - 1 .solve 



25 



N-l 

X r ik( n )Pkj( n ) = 



Vie [j,N-l] 



5:: 



j=0 



j-1 



k = 0 



for { p kj (n), Vke [j, N - 1] } 



30 (Equation 50) 

Note that the right hand sides of Equation (50) for all i 
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at each recursion step j do not contain any unknowns, 
i.e., {pj k (n)) there have already been found in previous 
•stages, Equation (45) and Equation (49) are just special 
cases of Equation (50), and {p k j(n), k=j,j+l,...,N-l) 
5 found in recursion step j form a column vector Pj(n) , which 
consists of the lower N-j elements of the j'th (0 < j < N- 
1) column of P(n) . The process of Equation (50) will take N 
divisions and 



10 



[2N 2 + 3N] + [2(N - 1 ) 2 + 4(N - 1 )] + [2(N - 2) 2 + 5(N - 2)] + ... 
+ [2(1) 2 + (N + 2)(1)] 

N 0 N _ N 

= ]T [2k" + (N + 3 - k)k] = ^k+(N + 3) 
k = l k = i k = i 

15 = N (N+1)(2N+1) + N (N+1)(N + 3)= | N (N+l)(N + 2) MACs 

(Equation 51) 

for steepest descent method, and 
N divisions and 

20 [2N 2 + 5N + 2] + [2(N - 1 ) 2 + 6(N - 1 ) + 2] + [2(N - 2) 2 + 7(N - 2) + 2] + ... 

+ [2(l) 2 + (N + 4)(l) + 2] 

N - N . N 

= £ [2k +(N + 5-k)k + 2] = k 2 + (N + 5) ]T k + 2N 

k = 1 k = 1 k = 1 



25 



= y (N+1)(2 N+l) + ^(N+l)(N + 5) + 2N= | N ( n2 + y N + y) MACs 



(Equation 52) 
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for conjugate gradient method. Note that in deriving Equa- 
tions (51) and (52) the following formulae are used 

]T k 2 = ~(N+ 1)(2N + 1) 
k = 1 

5 A N 

£ k = |(N + 1) 

k = 1 

(Equation 53) 

which can be easily proven by mathematical induction. 

Based on the above derivations, the SFAP method 
according to the third embodiment of the invention 
includes the following steps: 
Initialization : 

0 



10 



15 W(0) = 0 , X(0) = 0 , tj(0) = 0 , R(0) = 5I , s(0). = 0, P(0) = 



20 



25 



30 



(Equation 54) 

Updating the adaptive filter coefficients in sampling 
interval n including the steps shown in Equation 55 below. 
Please, note that designations used in Equation (55), are 
as follows: £(n) is defined in Equation (23) above, R(n) is 
the first column of R(n), R(n) is an N-l vector that con- 
sists of the N-l lower most elements of the N vector R(n), 
and fj(n) is an N-l vector that consists of the N-l upper 
most elements of the N vector rj(n) . Please, also note that 
any division operation in the 2nd expression of Equation 
(55) is not performed if the denominator is not greater 
than zero, in which case a zero is assigned to the quo- 
tient. 
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10 



15 



20 



25 



R(n) = R(n- l) + §(n)§ (n) 
-^(n-L)| T (n-L) 



P(n) = 



R '(n) 



Recursion with 
decrementing order 
to do matrix inversion 



j 



W(n) = W(n- U + ariN.^n- l)X(n-N) 



y(n) = W T (n)X(n) + ag T (n- l)R(n) 



e(n) = d(n)-y(n) 



e(n) 
(1 -a)e(n- 1) 



e(n) = 
e(n) = P(n)e(n) 



H(n) = 



0 

fj(n-l) 



+ e(n) 



Total 



MAC 



2N 



Division 



^N(N+ l)(N + 2) 
6 

(Steepest descent) 
or 



N 2 + 21 N + 28 



5 -n(: 

6 1 5 5 
(Conjugate gradient) 



L + N- 1 



N - 1 



2L + -N 3 + ^N 2 + — N-2 

6 2 3 
(Steepest descent) 

or 

2L + ?N 3 + -N 2 + — N-2 

6 2 3 
(Conjugate gradient) 



N 



N 



(Equation 55) 

An adaptive filter 300 according to a third embodiment of 
30 the invention, shown in Figure 5, is similar to that of 
Fig. 2 with like elements being designated by same refer- 
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ence numerals incremented by 200. The filter 300 also dif- 
fers from the filter 100 by the following features: the 
normilized step size may have any value from 0 to 1.0, the 
calculator 308 now has more extended structure for consec- 
5 utively determining columns of the inverse auto-correla- 
tion matrix in accordance with the steepest descent 
technique, and an e(n) calculator 320 is added. 

The P(n) calculator 308, n'ow being a matrix calcula- 
tor, operates in accordance with the flow-chart 400 shown 
10 in Figure 6. Upon start up for the sampling interval n 
(block 401) , the routine 402 sets an initial value to 
index j (block 404) which is submitted together with the 
• auto-correlation matrix R(n) (block 406) to a projection 
coefficient column calculator (block 408) . The calculator 
15 provides a steepest descent iteration in accordance with 
Equation (50) for the current value of index j, thus updat- 
ing the corresponding column of projection coefficients 
from the previous sampling interval (block 408) . The 
updated column of the projected coefficients is sent to a 
20 storage means (routine 410, block 412) to be stored until 
the other columns of P(n) are calculated. Until the index j 
is equal to N-l (block 416), its value is incremented by 1, 
i.e. made equal to j+1 (block 418), and the steepest 
descent iteration is repeated (block 408) to determine the 
25 next column of P(n). By performing N corresponding steep- 
est descent iterations f or j = 0, 1, ... N-l , all columns of the 
inverse auto-correlation matrix are thus determined and 
assembled into P(n) in an assembling means (block 414). A 
command/ signal (block 420) then notifies about the end of 
30 the sampling interval n and the beginning of the next sam- 
pling interval n+1 where the steps of the routine 4 00 are 
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repeated. In Figure 6, thick lines represent the propaga- 
tion of a matrix or vector signal, i.e., with more than 
one component, and the use of a thin line stands for a 
control propagation. 
5 In modification to this embodiment, the steepest 

descent calculator 308 may be replaced with the conjugate 
calculator. The corresponding structure is illustrated by 
a flow-chart 500 shown in Figure 7 where the blocks simi- 
lar to that ones of Figure 6 are designated by same refer- 
10 ence numerals incremented by 100. It operates in a manner 
described above with regard to Figure 6. 

A method of adaptive filtering according to a fourth 
embodiment of the present invention also provides adaptive 
filtering when the normalized step size has any value from 
15 0 to 1. It updates the adaptive filter coefficients by 
iteratively solving a number of systems linear equations 
which avoid an explicit matrix inversion performed in the 
third embodiment of the invention. The details are 
described below. 
20 The second equation from the set of Equations (6), 

which is reproduced for convenience in Equation (56) 
below, is equivalent to 

R(n)£(n) = £(n) 

(Equation 56) 

25 It is possible to obtain £(n), required for updating the 
adaptive filter coefficients, directly from the set of 
linear Equations (56) , which are solved again by one of 
the descending iterative methods. 

As a way of example, we will use a conjugate gradient 

30 method and perform N conjugate gradient iterations so that 
an exact solution, not an iterated one, is reached. It is 
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ensured by ■ the fact that the conjugate gradient method is 
guaranteed to reach the solution in no more than N itera- 
tions, with N being the order of the problem, see Equation 
(55) . It is convenient to start with £(n)=0 before itera- 
5 tions begin at each sampling interval n to save some com- 
putation time . 

Accordingly, the SFAP method of the fourth embodiment 
of the invention includes the following steps: 



10 



15 



20 



Initialization: 


MAC 


Division 


e t (0) = 0 , s(0) = 0 , r srs (0) = 0 , b(0) 


= 0 






In sampling internal n, repeat the following 
equations N times, i.e., for k = 0, 1, N-l: 






g = R(n)8 t (k)-e(n) 


' (N-l)xN 2 




J = r srs (k)g T b(k) 


(N — 1) x (N + 1 ) 




-— §(-k-+ 1-) = Y§(k)— g 










■ (N - 1 ) x N " 




b(k+ 1) = R(n)s(k+ 1) 


NxN 2 




r srs (k+ 1) - 

s(k+l)b(k+l) 


NxN 


Nx 1 


P = -r srs (k+l)g T s(k+l) 


Nx(N + 1) 




e t (k + 1) = e t (k) + P§(k+ 1) 


NxN 




Output: 






g(n) = e t (N) 








Total 


2N 3 +4N 2 - 1 


N 



30 



(Equation 57) 
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The steps of the adaptive filtering methods according 
to the fourth embodiment are presented in more detail 
below: 

Initialization : 

0 



W(0) = 0 , X(0) = 0 , r|(0) = 0 , R(0) = 81 , e(0) = 0 , P(0) = 



10 



15 



20 



25 



Processing in sampling interval n: 

R(n) = R(n- l) + |(n)§ T (n) 
-|(n-L)| T (n-L) 

W(n) = W(n- 1) + an N _,(n- l)X(n-N) 



y(n) = W T (n)X(n) + afj T (n- l)R(n) 



e(n) = d(n)-y(n) 



e(n) = 



e(n) 
(l-a)e(n- 1)J 



N conjugate 

gradient 

iterations 



( \J 
1 Solve R(n)e(n) = e(n) for e(n) 

S I 



r|(n) = 



0. 

Ti(n-l) 



+ e(n) 



Total 



(Equation 58) 



MAC 



2N 



L + N- 1 



N- 1 



2N 3 + 4N 2 - 1 



Division 



2L + 2N 3 + 4N 2 + 4N-3 



N 



N 



30 



(Ecruation 59) 

where the designations are similar to that presented with 
regard to the first, second and third embodiments 
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described above. Note that any division operation in Equa- 
tion (56) is not performed if the denominator is not 
greater than zero, in which case a zero is assigned to the 
quotient . 

5 An adaptive filter 600 according to a fourth embodi- 

ment of the invention is shown in detail in Figure 8. It 
includes a filter 602 characterized by adaptive filter- 
coefficients W(n), and means 604 for updating the coeffi- 
cients, the means being set with a normalized step size a 
10 having any value in a range from 0 to 1.0. The filter 602 
is a finite impulse response (FIR) filter which receives a 
reference input signal x(n) and an auxiliary signal f(n) 
used for updating the coefficients, and generates a provi- 
sional echo estimate signal PR(n) . The updating means 604 
15 includes a correlator 606 for recursively determining an 
auto-correlation signal presented in the form of auto-cor- 
relation matrix coefficients R(n) based on the reference 
input signal x(n), ari £(n) calculator 608 and an e(n) calcula- 
tor 62 0 for corresponding calculation of vectors £(n) and 
20 e(n). The calculator 608 defines £(n) by using an iterative 
conjugate gradient method having an inherent stability of 
operation as illustrated in detail above. The projection 
coefficients are used within updating means 604 for gener- 
ation the auxiliary filter adaptation signal f(n) and an 
25 echo estimate correction signal EC(n) . The latter is used 
together with the provisional echo estimate PR(n) to pro- 
duce the echo estimate signal y(n) . In Fig. 8 thick lines 
represent propagation of a matrix or vector signal, i.e., 
the signal with more than one component, and the use of a 
30 thin line stands for a scalar signal propagation. In Fig. 
8 a correlator 606 determines the autocorrelation matrix 
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R(n) in accordance with the first formula of Eq. (59) using 
the current and past x(n) samples. An n li(n) calculator" 610 
calculates n(n). based the last formula of Eq. (59) , and as 
shown in Fig. 8, n(n) is not used by the updating means 104. 
until the next sampling interval. The filter 602 produces 
the convolutional sum W A (n)X(n) . T| N _ 1 (n-l) is obtained from 
'Hn-iO) by putting the latter through a unit delay element 
611, providing a delay of one sampling interval, and fur- 
ther multiplied by the step size a in a Multiplier 613. 
The result is used for updating the adaptive filter coef- 
ficients (Eq. 59, second formula). n T (n-l) is dot-multiplied 
with part of R(n) by a Dot multiplier 612, and the result 
is further multiplied by a multiplier 614 with the step 
size a to form the correction term to be added to W T (n)X(n) 
by the summer 616 to form the filter output y(n) (Equation 
(59), third formula). Signals y(n) and e(n) are further sent 
to the e(n) calculator 62 0 to determine e(n) in accordance 
with a fourth and fifth formulae of Equation (59) , and the 
results are sent to the £(n) calculator 608 together with 
the auto-correlation matrix R(n) derived in the correlator 
606. The £(n) calculator 608 solves the sixth equation of 
Eq. (59) for £(n) by a conjugate gradient method, thus pro- 
viding sufficient data for updating the adaptive filter 
coefficients (Eq. 6, first formula) . 

The £(n) calculator 608, shown in detail in Figure 9, 
includes a one-step calculator 708a similar to the calcu- 
lator 2 08 of Fig. 4 and includes like elements which are 
referred to by the same reference numerals incremented by 
500 respectively (except for P(n-l) and Pfn) being replaced 
with £(n-l) and £(n) respectively) . Thick lines represent the 
propagation of a matrix or vector signal, i.e., with more 
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than one component, and the use of a thin line stands for 
a scalar signal propagation. At each sampling interval n, 
the calculator 708a performs N steps corresponding to k = 0, 
1,... N-l , each step being similar to the conjugate gradient 
5 iteration performed by the filter 208 of the second embod- 
iment of the invention. The calculator 608 additionally 
includes an output switch 754 which automatically opens- 
at the beginning of the sampling interval and closes at 
the end of N conjugate gradient iterations. 
10 Modifications described with regard to the first two 

embodiments are equally applicable to the third and fourth 
embodiments of the invention. 

Two " C M prototypes according to the third and fourth 
embodiments of the invention have been implemented in a 
15 floating point PC platform. They have demonstrated results 
completely consistent with the results of the first and 
second embodiments of the invention. 

Thus, an adaptive filter and a method providing a sta- 
bility of adaptive filtering based on feedback adjustment, 
20 are provided. 

Although the methods operate with real-valued num- 
bers, it does not prevent the invention from being 
extended to cases where introduction of complex numbers is 
necessary. 

25 Although the embodiments are illustrated within the 

context of echo cancellation, the results are also appli- 
cable to other adaptive filtering applications. 

Thus, it will be appreciated that, while specific 
embodiments of the invention are described in detail 

30 above, numerous variations, modifications and combinations 
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of these embodiments fall within the scope of the inven- 
tion as defined in the following claims. 
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WHAT IS CLAIMED IS : 

1. A method of adaptive filtering, comprising the 

steps of : 

(a) determining adaptive filter coefficients; 
5 (b) defining, a normalized step size; 

(c) updating the filter coefficients, comprising: 
determining auto-correlation matrix coefficients from 
a reference input signal, and 

solving at least one system of linear equations whose 
10 coefficients are the auto-correlation matrix coefficients, 
the system being solved by using a descending iterative 
method having an inherent stability of its operation, the 
results of the solution being used for updating the filter 
coefficients and the number of systems of linear equations 
15 to be solved being dependent on the normalized step size; 

(d) repeating the steps (b) and (c) required number 
of times . 

2. A method as defined in claim 1 wherein the step 

20 of determining auto-correlation matrix coefficients com- 
prises calculating the auto-correlation matrix coeffi- 
cients recursively. 

3. A method as defined in claim 1, wherein the 

25 step of defining a normalized step size comprises setting 
the normalized step size not equal to unity. 



30 



4. A method as defined in claim 1, wherein the 

step of defining a normalized step size comprises setting 
the normalized step size substantially less than unity. 
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5. A method as defined in claim 1, wherein the 

step of defining a normalized step size comprises setting 
the normalized step size less than about 0.7. 



5 6. A method as defined in claim 3, wherein the 

step of solving at least one system of linear equations 
comprises solving N systems of linear equations, with N 
being a projection order. 

10 7. A method as defined in claim 1, wherein the 

step of defining a normalized step size comprises setting 
the normalized step size close to unity. 

8. A method as defined in claim 1, wherein the 
15 step of defining a normalized step size comprises setting 

the normalized step size equal to unity. 

9. A method as defined in claim 1, wherein the 
step of defining a normalized step size comprises setting 

20 the normalized step size in a range from about 0.9 to 1.0. 

10. A method as defined in claim 1, wherein the 
step of defining a normalized step size comprises setting 
the normalized step size in a range from about 0.7 to 1.0. 

25 

11. A method as defined in claim 7, wherein the 
step of solving at least one system of linear equations 
comprises solving one system of linear equations only. 



30 12. A method as defined in claim 1, wherein the step 

of solving the system of linear equations by a descending 
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iterative method comprises solving the system by using a 
steepest descent method. 

13. a method as defined in claim 1, wherein the 
5 step of solving the system of linear equations by a 

descending iterative method comprises solving the system 
by using a conjugate gradient method. 

14. a method as defined in claim 1, wherein the 
10 step of solving the system of linear equations by a 

descending iterative method comprises solving the system 
by using a Newton ' s method . 

15. a method as defined in claim 1, wherein the 
15 step of solving the system of linear equations by a 

descending iterative method comprises solving the system 
by using a PARTAN method. 

16. a method as defined in claim 1, wherein the 

20 step of solving the system of linear equations by a 
descending iterative method comprises solving the system 
by using a quas i -Newton ' s method. 

17. A method as defined in claim 1, wherein the 
25 steps are performed by operating with real value numbers. 

18. A method as defined in claim 1, wherein the 
steps are performed by operating with complex value num- 
bers . 



30 



19. A method as defined in claim 1, the method 
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being used in an application selected from the group con- 
sisting of echo cancellation, noise cancellation, channel 
equalization and system identification. 

5 20. A method as defined in claim 1, wherein the 

step of solving at least one system of linear equations 
comprises determining projection coefficients, the projec- 
tion coefficients being the coefficients of an inverse 
auto-correlation matrix. 



10 



15 



21. A method as defined in claim 20, wherein the 
step of determining auto-correlation matrix coefficients 
comprises calculating the auto-correlation matrix coeffi- 
cients recursively . 

22. A method as defined in claim 20, wherein the 
step of defining a normalized step size comprises setting 
the normalized step size not equal to unity. 



20 23. A method as defined in claim 20, wherein the 

step of defining a normalized step size comprises setting 
the normalized step size substantially less than unity. 

24. A method as defined in claim 20, wherein the 
25 step of defining a normalized step size comprises setting 

the normalized step size less than about 0.7. 

25. A method as defined in claim 22; wherein the 
step of solving at least one system of linear equations 

30 comprises solving N systems of linear equations, with N 
being a projection order. 
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26. A method as defined in claim 25, wherein the 
step of solving N systems of linear equations comprises 
solving N systems of linear equations having decrementing 

5 orders . 

27. A method as defined in claim 20, wherein the 
step of defining a normalized step size comprises setting 
the normalized step size close to unity. 

10 

28. a method as defined in claim 20, wherein the 
step of defining a normalized step size comprises setting 
the normalized step size equal to unity. 

15 29. A method as defined in claim 20, wherein the 

step of defining a normalized step size comprises setting 
the normalized step size in a range from about 0.9 to 1.0. 

30. A method as defined in claim 20, wherein the 
20 step of defining a normalized step size comprises setting 

the normalized step size in a range from about 0.7 to 1.0. 

31. A method as defined in claim 27, wherein the 
step of solving at least one system of linear equations 

25 comprises solving one system of linear equations only. 

32. A method as defined in claim 31, wherein 

determining the projection coefficients comprises calcu- 
lating coefficients of a first column of the inverse auto- 
30 correlation matrix coefficients only. 
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33. A method as defined in claim 20, wherein the 
step of solving the system of linear equations by a 
descending iterative method comprises solving the system 
by using a steepest descent method. 

5 

34. A method as defined in claim 20, wherein the 
step of solving the system of linear equations by a 
descending iterative method comprises solving the system 
by using a conjugate gradient method. 
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35. A method as defined in claim 20, wherein the 

step of solving the system of linear equations by a 
descending iterative method comprises solving the, system 
by using a Newton's method. 

36. A method as defined in claim 20, wherein the 

step of solving the system of linear equations by a 
descending iterative method comprises solving the system 
by using a PART AN method. 

37 m a method as defined in claim 20, wherein the 

step of solving the system of linear equations by a 
descending iterative method comprises solving the system 
by using a quasi-Newton's method. 

33 A method as defined in claim 20, wherein the 

steps are performed by operating with real value numbers. 

39^ a method as defined in claim 20, wherein the 

30 steps are performed by operating with complex value num- 
bers . 
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40. A method as defined in claim 20, the method 
being used in an application selected from the group con- 
sisting of echo cancellation, noise cancellation, channel 

5 equalization and system identification. 

41. An adaptive filter, comprising: 

a filter characterized by adaptive filter coeffi- 
cients ; 

10 means for updating the filter coefficients, including 

means for setting a normalized step size, the updating 
means comprising: 

a correlator for determining auto-correlation matrix 
coefficients from a reference input signal, and 
15 a calculator for solving at least one system of lin- 

ear equations whose coefficients are the auto-correlation 
matrix coefficients, the system being solved by using a 
descending iterative method having an inherent stability 
of its operation, the results of the solution being used 
20 for updating the filter coefficients and the number of 
systems of linear equations to be solved being dependent 
on the normalized step size. 

42. The adaptive filter as defined in claim 41, 
25 wherein the correlator is a recursive correlator . 

43. The adaptive filter as defined in claim 41, 
wherein the normalized step size is not equal to unity. 



30 44. The adaptive filter as defined in claim 41, 

wherein the normalized step size is substantially less 
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than unity. 

45. The adaptive filter as defined in claim 41, 

wherein the normalized step size is less than about 0.7. 

5 

46. The adaptive filter as defined in claim 43, 
wherein the calculator includes means for solving N sys- 
tems of linear equations, with N being a projection order. 

10 47. The adaptive filter as defined in claim 41, 

wherein the normalized step size is close to unity. 

48. The adaptive filter as defined in claim 41, 
wherein the normalized step size is equal to unity. 

15 

49. The adaptive filter as defined in claim 41, 
wherein the normalized step size is within a range from 
about 0.9 to 1.0. 

20 50. The adaptive filter as defined in claim 41, 

wherein the normalized step size is within a range from 
about 0 . 7 to 1.0. 

51. The adaptive filter as defined in claim 47, 

25 wherein the calculator provides solution of one system of 
linear equations only. 

52. The adaptive filter as defined in claim 41, 
wherein the calculator is a calculator providing solution 
30 of the system of linear equations according to a steepest 
descent method . 
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53. The adaptive filter as defined in claim 41, 
wherein the calculator is a calculator providing solution 
of the system of linear equations according to a conju- 

5 gate gradient method . 

54. The adaptive filter as defined in claim 41,- 
wherein the calculator is a calculator is a calculator 
providing solution of the system of linear equations 

10 according to a Newton's method. 

55. The adaptive filter as defined in claim 41, 
wherein the calculator is a calculator providing solution 
of the system of linear equations according a PARTAN 

15 method. 

56. The adaptive filter as defined in claim 41, 

wherein the calculator is a calculator providing solution 
of the system of linear equations according to a quasi- 
20 Newton's method. 

57. The adaptive filter as defined in claim 41 

capable of operating with real value numbers. 



25 



58. The adaptive filter as defined in claim 41 

capable of operating with complex value numbers. 



59. The adaptive filter as defined in claim 41 for 

use in an application selected from the group consisting 
30 of echo cancellation, noise cancellation, channel equal- 
ization and system identification. 
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60. The adaptive filter as defined in claim 41, 
wherein the calculator further comprises means for deter- 
mining projection coefficients, the projection coeffi- 

5 cients being the coefficients of an inverse auto- 
correlation matrix. 

61. The adaptive filter as defined in claim 60, 
wherein the correlator is a recursive correlator. 

10 

62. The adaptive filter as defined in claim 60, 
wherein the normalized step size is not equal to unity. 

63.. The adaptive filter as defined in claim 60, 

15 wherein the normalized step size is substantially less 
than unity. 

64. The adaptive filter as defined in claim 60, 
wherein the normalized step size is less than about 0.7. 

20 

65. The adaptive filter as defined in claim 62, 
wherein the calculator is capable of solving N systems of 
linear equations, with N being a projection order. 

25 66. The adaptive filter as defined in claim 65, 

wherein the calculator is capable of solving N systems of 
linear equations having decrementing orders. 



67. The adaptive filter as defined in claim 60, 

30 wherein the normalized step size is close to unity. 
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68. The adaptive filter as defined in claim 60, 
wherein the normalized step size is equal to unity. 

69. The adaptive filter as defined in claim 60, 
5 wherein the normalized step size is within a range from 

about 0.9 to 1.0.- 

70. The adaptive filter as defined in claim 60, 
wherein the normalized step size is within a range from 

10 about 0.7 to 1.0. 

71. The adaptive filter as defined in claim 67, 
wherein the calculator is suitable for solving one system 
of linear equations only. 

15 

72. The adaptive filter as defined in claim 71, 
wherein the means for determining projection coefficients 
provides calculation of coefficients of a first column of 
the inverse auto-correlation matrix coefficients only. 

20 

73. The adaptive filter as defined in claim 60, 
wherein the calculator is a calculator providing solution 
of the system of linear equations according to a steepest 
descent method. 

25 

74. The adaptive filter as defined in claim 60, 
wherein the calculator is a calculator providing solution 
of the system of linear equations according to a conjugate 
gradient method. 

30 

75. The adaptive filter as defined in claim 60, 
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wherein the calculator is a calculator providing solution 
of the system of linear equations according to a Newton's 
method . 

5 76. The adaptive filter as defined in claim 60, 

wherein the calculator is a calculator providing solution 
of the system of linear equations according to a PARTAN 
method . 

10 77 . The adaptive filter as defined in claim 60, 

wherein the calculator is a calculator providing solution 
of the system according to a quasi -Newton ' s method. 

78. The adaptive filter as defined in claim 60 
15 capable of operating with real value numbers. 

79. The adaptive filter as defined in claim 60 
capable of operating with complex value numbers. 

20 80. The adaptive filter as defined in claim 60 for 

use in an application selected from the group consisting 
of echo cancellation, noise cancellation, channel equal- 
ization and system identification. 

25 
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